function mdl = estimate(mdl, X, w, copts, cachelevel)
%ESTIMATE Estimates the Gaussian model(s) from data
%
% [ Syntax ]
%   - mdl = estimate(mdl, X, [])
%   - mdl = estimate(mdl, X, w)
%   - mdl = estimate(mdl, X, w, options)
%   - mdl = estimate(mdl, X, w, copts, cachelevel)
%
% [ Arguments ]
%   - mdl:          The Gaussian model(s)
%   - X:            The sample data
%   - w:            The weights of the samples
%   - options:      The string giving the options in compact representation.
%   - cachelevel:   The level of pre-computation cache
%
% [ Description ]
%   - mdl = estimate(mdl, X, []) estimates a Gaussian model based 
%     on the samples given in X. 
%
%   - mdl = estimate(mdl, X, w) estimates Gaussian model(s) based
%     on the samples given in X. 
%
%     If w is a row vector, then a single model is estimated. If w is
%     a matrix of size K x n, then K models will be estimated, and a
%     object array will be returned.
%
%   - mdl = estimate(mdl, X, w, options) estimates Gaussian model(s)
%     with further options. 
%
%     The options are specified with compact string represenation. 
%     Here are the available options:
%       \{:
%           - 'sharecov':    all models share the same covariance
%           - 'fixmean':     the mean vectors of the models are fixed
%                            to those in the input mdl
%           - 'fixcov':      the covariance matrices of the models are
%                            fixed to those in the input mdl
%       \:}
%     By default, these options are set off. To open them, you may 
%     specify them in the options. For example, if the options are
%     'sharecov', then the covariance-sharing is turned on; if the
%     options are 'sharecov|fixmean' then both sharecov and fixmean
%     are turned on. (Options are separated by |).
%
%   - mdl = estimate(mdl, X, w, copts, cachelevel) estimate the models
%     and set their cache to specified level. 
%
%     By default, the cachelevel is set to 3.
%
% [ Remarks ]
%   - Generally, the input mdl only serves as an agent such that a
%     proper estimate method is called. So it is typically an empty
%     model.
%
%   - If either fixmean or fixcov is turned on, then in the input,
%     mdl should not be empty, such that the mean vectors or covariance
%     matrices are ready there. 
%
% [ History ]
%   - Created by Dahua Lin, on Dec 23, 2007
%

%% parse and verify input arguments

assert(isnumeric(X) && ndims(X) == 2, ...
    'sltoolbox:slcfullgauss:estimate:invalidarg', ...
    'X should be a numeric matrix.');

if nargin < 3
    w = [];
end

if isempty(w)
    K = 1;
else
    assert(isnumeric(w) && ndims(w) == 2 && size(w,2) == size(X,2), ...
        'sltoolbox:slcfullgauss:estimate:invalidarg', ...
        'w should be a numeric matrix with n columns.');
    K = size(w, 1);
end

opts = struct( ...
    'sharecov', false, 'fixmean', false, 'fixcov', false);

if nargin >= 4
    opts = slcompactopts(opts, copts);
end

if opts.fixmean || opts.fixcov
    
    assert(K == numel(mdl), ...
        'sltoolbox:slcfullgauss:estimate:invalidarg', ...
        'The number of models in input mdl should be equal to K when fixmean or fixcov is on,');    
    
    for i = 1 : K
        assert(~isemptymdl(mdl(i)), ...
            'sltoolbox:slcfullgauss:estimate:emptymdl', ...
            'All models should not be empty when fixmean or fixcov is on.');
    end
    
    assert(~(opts.fixmean && opts.fixcov), ...
        'sltoolbox:slcfullgauss:estimate:invalidarg', ...
        'fixmean and fixcov cannot be both turned on');
end


if nargin < 5
    cachelevel = 3;
else
    assert(isnumeric(cachelevel) && isscalar(cachelevel) && cachelevel >= 0, ...
        'sltoolbox:slcfullgauss:estimate:invalidarg', ...
        'cachelevel should be a non-negative scalar.');
end



%% main

if ~opts.fixmean && ~opts.fixcov  
    
    % estimate a brand-new model
    
    % estimate mean    
    mu = slmean(X, w);
    
    % estimate    
    if ~opts.sharecov || K == 1
        sigma = slcov(X, w, mu);
    else
        sigma = {slpoolcov(X, w, mu)};
    end
    
    % construct model
    mdl = slcfullgauss(mu, sigma, cachelevel);
    
else
    
    % estimate based on input model
   
    if ~opts.fixcov     % re-estimate the covariance with mean fixed
       
        if ~opts.sharecov
            sigma = slcov(X, w, [mdl.mu]);
            for i = 1 : K
                mdl(i).sigma = sigma(:,:,i);
            end
        else
            sigma = slpoolcov(X, w, [mdl.mu]);
            for i = 1 : K
                mdl(i).sigma = sigma;
            end
        end        
        
    elseif ~opts.fixmean  % re-estimate the mean with cov fixed
        
        mu = slmean(X, w);        
        for i = 1 : K
            mdl(i).mu = mu(:,i);
        end
        
    end    
    
    % re-cache
    mdl = cache(mdl, 0);
    
    if cachelevel > 0
        mdl = cache(mdl, cachelevel);
    end    
end


